Cumulants of heat transfer in nonlinear quantum systems 
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We consider thermal conduction across a general nonlinear phononic junction. Based on two- 
time observation protocol and the field theoretical method, the cumulants of the heat transfer in 
both transient and steady-state regimes are studied on an equal footing, and a practical formula 
for cumulant generating function of heat transfer is obtained. To demonstrate the power of the 
general formalism, we study anharmonic effects on fluctuation of steady-state heat transfer across 
a typical molecular junction with a quartic nonlinear on-site pinning potential. In addition, an 
explicit nonlinear modification to the cumulant generating function exact up to the first order of 
nonlinear strength is given. Gallavotti-Cohen fluctuation symmetry is also verified. As a numerical 
illustration to first few cumulants of heat transfer, a self-consistent procedure is introduced, which 
works well for strong nonlinearity. 
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The physics of nonequilibrium many-body systems is 
one of the most rapidly expanding areas which attracts 
much attention recently. With the development of the 
modern nanoscale technology, a full understanding of the 
general features of thermal conduction is needed. And it 
is well-known that the noise [l[ generated by nanodevices 
contains valuable information on microscopic transport 
processes not available from current, which is simply re- 
lated to the lowest order moment of the full distribution 
function for heat transfer. To this end, the full counting 
statistics (FCS) regarding the distribution of transferred 
quantity, such as heat in phononic case, has to be deter- 
mined. 

The extensive study of the FCS started from the field 
of electronic transport pioneered by Lcvitov and Leso- 
vik Q, while much less attention is given to phonon 
transport. Saito and Dhar are the first ones to borrow 
this concept to thermal transport However, both of 
them, including many subsequent works, are mainly re- 
stricted to noninteracting systems 0, [B| , although some 
works have already been devoted to the analysis of fluctu- 
ation considering the effect of nonlinearity in the classical 
limit through Langcvin simulations, or approximately in 
a restricted electronic transport case, such as the FCS in 
molecular junctions with electron- phonon interaction Q . 

In recent years, phononics, i.e., the counterpart tech- 
nology of electronics, has great interest to both theorists 
and experimentalists, and presents an unforeseen wealth 
of applications 0, H| . And the nonlinearity, such as the 
phonon-phonon interaction, has been found of special im- 
portance in the construction of phononic devices In 
addition, recently it has been noted that the nonlinearity 
of interaction is crucial to the manifestation of geomet- 
ric heat flux [§]. Therefore, it is desirable to establish a 
systematic and practical formalism to properly deal with 



cumulants of heat transfer in the presence of nonlinearity. 

In this work, we shall study the FCS for heat trans- 
fer flowing across a quantum junction in the presence 
of general phonon-phonon interactions. Then, based on 
two-time observation protocol (lp| . we construct a con- 
cise and rigorous cumulant generating function (CGF) 
expression for the heat transfer in this general situation, 
in which both the transient and steady-state case are 
dealt with on equal footing. Furthermore, as an illus- 
tration of this general formalism, a single site junction 
with a quartic nonlinear on-site pinning potential is in- 
troduced and the corresponding CGF exact up to first 
order of the nonlinear interaction strength is given. Fi- 
nally, a self-consistent scheme is employed to numerically 
illustrate the nonlinear effects on first three cumulants of 
heat transfer. 

Formalism — We consider the lead-junction-lead 
model initially prepared in a product state pi„i = 



n a =L,c,.R Tr '( e -g g -H g ) ■ It can be imagined that left lead 
(L), center junction (C), and right lead (R) in this model 
were in contact with three different heat baths at the in- 
verse temperatures /3l — {ksTh) 1 , fic — (fcsTc) -1 and 
Pr = {ksTji)" 1 , respectively, for time t < t . At time 
t = to, all the heat baths are removed, and couplings of 
the center junction with the leads rife = u I ]V LC uc and 



H 



CR 



•c 



V iir and the nonlinear term rl n appearing 



in the center junction are switched on abruptly. Now the 
total Hamiltonian is given by 



Hi 



Hl + Hc 



Hlc + H 



CR 



(1) 



where rl a = ^p^p Q +^u^K a u a , a = L,C,R, repre- 
sents coupled harmonic oscillators, and 
p a are column vectors of transformed coordinates and 
corresponding conjugate momenta in region a. The su- 
perscript T stands for matrix transpose. The formalism 



2 



developed here does not require to specify the explicit 
form of H n . 

Based on two-time observation protocol, and employ- 
ing the consistent quantum framework 11] , we can obtain 
the generating function (GF) for heat transfer during the 
time tu — to to be 



a . b 



(2) 



where, Pr (P t " P t b M ) stands for the joint probability for 
the quantum history P t " P t b M , that is, the result of 
the measurement at time to of energy of the left lead 
associated with the operator Hl is the eigenvalue a of 
Hl, then the measurement at time tu yields the eigen- 
value b of Hl', M-£/2 (tM,to) means evolution operator 
associated with the modified total Hamiltonian H t ^ t = 
e <(-f/2)'Hi^ tote -i(-e/2)«i j S i m ii ar iy for% 2 (*o,*m) and 
(. . .) denotes the ensemble average over the initial state 

Pint • 

Transforming to the interaction picture with respect 
to the free Hamiltonian h = Hl + He + Hr, the GF 
can be written as Z (f) = (T T e~T, fc ^T ? (r)^ ^ wne re T T 
is a r-ordering operator arranging operators with the 
earliest r on the contour C (from to to £m and back 



to to) to the right, and 7^ (r) 



(r) V LC u c (r) + 



«P (t) V ur (r) + "Hn (t) (a caret is put above oper- 
ators to denote their r dependence with respect to the 
free Hamiltonian h such as uc ( T ) = es McT uce"^ Wcr ), 
and u L (r) = ul {frx T + t) with x T = — £/2 (£/2) when 
t = t + (i~ ) is on the upper (lower) branch of the contour 
C. 

Employing the field theoretical/diagrammatic method 
[l2^ , it can be shown that the derivatives of the connected 
vacuum diagrams are closely related to the Green's func- 
tions. Specifically, we could obtain 



dlnZ _1 
in terms of 



drdr'Ti 



c 



Gcc (r,r') 



d± L (t',t) 

d(iO 



(3) 



Gcc (t~i, t 2 ) 



(t t u c .j (n) u c .k (r 2 ) e-i ^ I 



(4) 



and £ L = V CL ~g L V LC with g L (n,T 2 ) jk = 
-k(TrUl d (Ti)ul tk (r 2 )) and y CL = (y ic ) T The 
tilde on the Green's functions emphasizes the fact that 
they are counting field ^-dependent. And if needed, the 
proper normalization for the CGF, i.e., lniJ(£), can be 
determined by the constraint In Z (0) = 0. The real-time 
version of Eq. which is essentially a generalized 



Meir-Wingreen formula for current, can be also shown 
algebraically using a similar method as Ref. [l3| based 
on nonequilibrium Fcynman-Hcllmann theorem. 

In order to evaluate In Z, a closed equation for 
Gcc { T i-, T 2) is needed, which could be deduced by a 
transformation from O (r) to O 1 (r) taking t£ as a refer- 
ence time such as u c (r) = V (tj , t) uc (t) V (r, tg") with 



W LC u c +i^V CR u R ] 



(5) 



and V(to~,r) = V (r, t$) , where the subscript 
C [t, £q ] denotes the path along the contour C from £q 
to t. The present transformation defined on the con- 
tour C is necessary, since the coupling Hamiltonian with 
counting field is different on the upper and lower branch 
respectively. Then according to Eq. @ , Gcc { T i, T 2) can 
be rewritten as 



Gcc (n, t 2 ) 



jk 



-Tr 



pL^,(r 1 )^(r 2 )e-t/ c ^(r) 



(6) 



in terms of p ini = p ini A/ Z and Z n = Z/Z , where 

z Q = (r T e-TJcH^L T v LC uc+u T c v CR u R )\j is the GF 

when H n = 0, A = V (*„ , tjif ) V (t M , ) , and Tr{p( m ) = 
1. Now the benefit of the transformation to the sec- 
ond interaction picture defined on the contour is clear, 
that is, the extra term A always appearing at the left- 
most position after contour-ordering can be taken out 
of T T to combine with p.i n i/Zo to yield pf ni . Still in 
this transformed picture / the Wick theorem is valid, 
and thus the Dyson equation for Gcc is obtained as 
Gcc — G° cc + G cc Y> n Gcc, a symbolic notation of 

Gcc (ti,t 2 ) = G cc (ti,t 2 ) 

+ / drdr'G cc (n,T)t n (T,r')Gcc (r' ,r 2 ) (7) 
Jc 



in terms of 



CC 



Pini T TU c (n) (r 2 ) 
9c + 9c \pL + ^r^J G cc 



(8) 
(9) 



and the nonlinear self energy S n constructed by Gq C 
solely due to the nonlinear Hamiltonian H n , where 
5jfl( T ii r 2) is the right- lead version of the ordinary 
contour-order self energy £„ = V Cv g v V vC , v = L,R, in 
which g a (T 1 ,T 2 ) jk = -\ (T T u aJ (n) u a>k (r 2 )) for a 
!..('. II are the uncoupled contour-order Green's func- 
tions. After introducing an auxiliary equation Gq C = 
9c + 9c + Si?) G^ c , and combining it with Eqs. ([7]) 
and (|9|), the closed Dyson equation for Gcc (ti, t 2 ) could 
be obtained as 



cc — (jc 



G 



cc ( 



+ E n ) G. 



'CC, 



(10) 
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where the shifted self energy = 11 l — Hl, which first 
appear in Ref. [f|, accounts for the FCS in ballistic sys- 
tems. 

From now on, for notational simplicity, all the sub- 
scripts GG of the Green's functions will be suppressed 
and the superscript in both G^ c and G° cc will be re- 
expressed MS i\ subscript. 

The above formalism is valid in both transient and 
steady case. Proceeding to study cumulants of steady- 
state heat transfer explicitly, one simply set to ~~ ► 
— oo and tM +00 simultaneously, and technically 
assume that real-time versions of G(t\,T2) are time- 
translationally invariant. Then going to the Fourier 



space, Eq. ((3|) for 

as 



am 



in steady state could be rewritten 



dlnZ 

Wo 



(tM — to) 



'dco 
n 27T 



hiu Tr 



G<H>e 



(11) 



after taking into account G > 



-uS\ = G < [uj] t and 
[uj\ . In the Fourier space, due to Eq. (JTUJ) 
exact result for G [u>] could be yielded as 



G[w] = (G M 1 - E A [w] - E„ [w] 



(12) 



when keeping in mind the convention that the contour- 
order Green's function such as G(ti,t 2 ) in frequency 
space is written as 



G| 



G* [u] G< [cj] 
-G>M -&[u] 



(13) 



For nonequilibrium Green's function (NEGF) notations 
and relations among Green's functions Go, we refer to 
Ref. 0. 

Application — Now we apply the general formalism 
developed above to study a monatomic molecule with 
a quartic nonlinear on-site pinning potential, that is, 
Hn = jXuq in Eq. ([T}. In this case, nonlinear contour- 
order self energy exact up to first order in nonlinear 
strength is E n (r, r') = 3ih\Go (r, t') S (r, r'), where the 
generalized S- function S (r, r') is the counterpart of the 
ordinary Dirac delta function on the contour G, see, for 
example, Ref. [lij]. Thus the corresponding frequency- 
space nonlinear self energy is 



S n [lj] = 3ih\ 



G* (0) 
G o (0). 



(14) 



Consequently, exact up to first order in nonlinear 
strength the CGF for the molecular junction could be 
given as 

1 <91nZ(£) 



(t M - h) d(i£) 
x GoHGo (0)-G§[w]G§ (0) 



[°° du; f d\nD[u} 



d 



d(i£) D[u 



(15) 



with 

D[ui] = dct [/-Go [ui]E A [uj] 

= 1 -T[w][{e«»» - l) h ( 1 + fn) + hu - 1) f R ( 1 + h) 

(16) 

and Gf{Q) = where T[u] = 

Tr (GqTjiGqTl) is the transmission coefficient in the bal- 
listic system, and /„ = [exp (j3 v hhj) — 1] , v = L, R 
is the Bosc-Einstein distribution function for phonons. 
Here G r = G - G^ and Gg = G^ - G& are retarded 
and advanced Green's functions, respectively. Also T v = 
i — E£] , v = L, R, related to the spectral density of 
the baths, are expressed by retarded and advanced self 
energies similarly defined as Green's functions. Eq. ([15]) 
satisfies Gallavotti-Cohen symmetry for the derivatives 
[l5| because D[cj] remains invariant under the transfor- 
mation £ — > — £ + — /?l) while <9D[u;]/<9(i£) changes 
sign. For A = an integration over £ can be performed 
trivially and the ballistic result can be recovered @, 0] . 

One could easily use this CGF in Eq. (Ti~5l) to evaluate 
cumulants. The steady current out of the left lead is 
closely related to the first cumulant so that 



JSS = 



d (dlnZ(£) 



dt M \ d(i£) 
dui 



e=o y 



47T 



kwT[u>] (1 + A[w])(/ L -/ B ) 



(17) 



where A[u] = 3^AG (0) (Gg + G r ) [uj] is the first-order 
nonlinear correction to the transmission coefficient. 

The fluctuation for steady-state heat transfer in the 
molecular junction is obtained by taking the second 
derivative with respect to i£, and then setting £ = 0: 



{{Q 2 )) 

{tM — to) 



^{(M 2 r 2 M(i + 2AM)(/ i -/ fl ) 



3/TAcj 



G t [u]SG t -G t [u]SGo T[u}(f L -f R ) 



(hxjf T[w](l + AH)(/i+/ s +2/ L / fl )}, (18) 



where, 



SGI 



d&o* (0) 



du) 
2^ 



huT[u] (/ £ -/ji)G* '*M 
(19) 



Higher-order cumulants can be also systematically given 
by corresponding higher-order derivatives. 

In the following, we will give a numerical illustra- 
tion to the first few cumulants for heat transfer in this 
molecular junction using a self-consistent procedure [lij ]. 
which means that the nonlinear contour-order self en- 
ergy is taken as E„(t,t') = 3ih\G (t, r') S (r, r'). We 
sclf-consistently calculate G [uj] , and higher deriva- 

tives, based on Eq. (JT2J) , then the first few cumulants are 
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X [eV/(amu 2 A 4 )] X [eV/(amu 2 A 4 )] 

FIG. 1. (Color online) The steady-state cumulants 
((Q n )) I '{pM — to) for n = 1,2,3 with nonlinear coupling 
strength A for k = 1 eV/(uA 2 ), k = 0.1 k, K c = 1.1 k, and 
V-°Q = Vo C i R = -0-25 k. The solid (dotted) line shows the 
self-consistent (first-order in A) calculation for the cumulants. 
The temperatures of the left and right lead are 660 K and 410 
K, respectively. 



obtained by the corresponding derivative of lnZ(£) in 
Eq. (fTT|) with respect to i£ and finally taking the value 
£ = 0. In order to obtain m-th order cumulant one needs 
to solve (m — l)-th derivative of G iteratively. For exam- 
ple, the iterative equation for computing the first deriva- 
tive of G is given as 



dGU 



f dt L [Lo] 



dt n [u] 



G[u 



(20) 

where G[u>] = G[cj]\^_ q . The equation transforms to 

a linear equation for dG[uj]/d(i£) by using the modi- 
fied £„ [&] , which is then solved numerically by iteration. 
Equations for the higher derivatives can be similarly ob- 
tained. Very recently, it is shown that such self-consistent 
calculation gives extremely accurate results for the cur- 
rent in the case of a single site model as compared with 
quantum master equation approach [l7j |. thus we believe 
that it also leads to excellent predictions for the higher 
order cumulants. Rubin baths are used for numerical cal- 
culation, that is, K a , a = L, R in Eq. ((T|) are both the 
semi-infinite tridiagonal spring constant matrices consist- 
ing of 2k + ko along the diagonal and — k along the two 
off-diagonals. And only the nearest interaction V_[ and 
Vqi 1 between the molecule and the two baths are con- 
sidered and He = \"Pc o + \^C^c o- 

Figure 1 shows the plot for the first three cumulants 
with nonlinear strength A. The effect of nonlinearity is 
to reduce the current as well as higher order fluctua- 
tions. For weak nonlinearity the results can be compared 
with the first-order perturbation ones, presented as dot- 
ted lines. The deviation for high A is understandable. 
For certain parameters (not shown) the third cumulant 
can change sign from positive to negative. The fact that 
third and higher order cumulants are small but nonzero 



implies that the distribution for transferred energy is not 
Gaussian. Similar effect was also observed for a classical 
system 0]. 

Summary — A formally rigorous formalism dealing 
with cumulants of heat transfer in nonlinear quantum 
thermal transport is established based on field theoret- 
ical and NEGF methods. The CGF for the heat trans- 
fer in both transient and steady-state regimes is stud- 
ied on an equal footing and useful formulas for the CGF 
are obtained. A new feature of this formalism is that 
counting-field dependent full Green's function can be ex- 
pressed solely through the nonlinear term fi^ (t) with an 
interaction picture transformation defined on a contour. 
Clearly, this method can be readily generalized to handle 
other nonlinear systems, such as electron-phonon interac- 
tion and Joule heating problems. Up to the first order in 
the nonlinear strength for the single-site quartic model, 
the CGF for steady-state heat transfer is obtained and 
explicit results for the steady current and fluctuation for 
steady-state heat transfer are given. Moreover, a self- 
consistent procedure is introduced to give a numerical 
illustration for first few order cumulants. It is found that 
nonlinearity generally suppresses fluctuations. 
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